The phase diagram of 2D polar condensates in a magnetic field 
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Spin one condensates in the polar (antiferromagnetic) phase in two dimensions are shown to undergo a tran- 
sition of the Ising type, in addition to the expected Kosterlitz-Thouless (KT) transition of half vortices, due to 
the quadratic Zeeman effect. We establish the phase diagram in terms of temperature and the strength of the 
Zeeman effect using Monte Carlo simulations. When the Zeeman effect is sufficiently strong the Ising and KT 
transitions meet. For very strong Zeeman field the remaining transition is of the familiar integer KT type. 
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Ultracold atomic gases represent a new frontier in quan- 
tum magnetism, where optical traps allow for the possibility 
of spontaneous ordering of the atoms' spin. In particular, sys- 
tems of spinful bosons allow for the possibility of studying 
the interplay of Bose condensation (or superfluidity) and mag- 
netic ordering in spinor condensates. In this Letter we discuss 
one of the simplest such systems: the polar (or antiferromag- 
netic) spin-1 condensate in a magnetic field, realized in a gas 
of 23 Na atoms (T|. We show that there are two types of de- 
fects: vortices and domain walls, or strings. Some years ago a 
number of authors lHJfil studied the interesting phase diagram 
that results from the competition between vortex interactions 
and string tension in a simple statistical mechanical model. 
Our goal here is to show that the same physics can arise in a 
atomic gas from quite different microscopic origins. 

Consider a dilute gas of spin-1 bosons, described by a 
spinor 0(r). The character of the magnetic order that develops 
at low temperature depends upon the interatomic interactions, 
described by two kinds of quartic terms (5]|6l 
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corresponding to a density-density and spin-spin interaction 
respectively. Here S are the spin-1 angular momentum ma- 
trices. It is convenient to work in terms of Cartesian compo- 
nents, where (Si)j k — —itijk- Writing = a + ih, where 
a and b are real vectors, the second term in Eq. ([T]) becomes 
2c 2 (a x by. We see that for c 2 < the energy is minimized 
at fixed density for a, b perpendicular and equal, while for 
the case c 2 > that will be our primary concern, a and b are 
parallel. Conventionally these two possibilities are called the 
ferromagnetic and polar states, respectively. 

In the mean-field description sketched above, a non-zero 
expectation value (0), or equivalently off-diagonal long- 
range order in the density matrix p a &(r, r') = (<^(r)</> b (r')), 
simultaneously describes Bose condensation and the breaking 
of rotational symmetry. It is natural to ask whether these two 
phenomena necessarily go hand-in-hand, and if not, which oc- 
curs first as the temperature is lowered. We will show that 
in the two dimensional polar system (c 2 > 0, appropriate 
to 23 Na), such an intermediate phase can arise, possessing 
quasi-long-range order in the singlet pair amplitude <f> ■ <p - 
a pair superfiuid (PS) - while the spin remains disordered. 
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FIG. 1: (Color online) The phase diagram for a two dimensional po- 
lar condensate with quadratic Zeeman effect. The dashed red line 
and solid blue line mark the KT transition and Ising transition re- 
spectively. In the normal phase (N) half vortices of equal or oppo- 
site charge are joined by domain walls (blue lines). In the superfiuid 
phase (S) vortices of opposite charge are bound (red ellipse) and there 
are no domain walls. In the pair superfiuid phase (PS), q/t < 4, vor- 
tices are bound but domain walls remain. 



Spin ordering occurs in a second transition at a lower tem- 
perature. This is possible when there is a small anisotropy 
originating from the quadratic Zeeman effect, permitting an 
Ising transition where none would be allowed at zero field by 
the Mermin-Wagner theorem. At larger anisotropies the PS 
phase vanishes. The resulting phase diagram, shown in Fig. 
[T] for our Monte Carlo simulations of a particular model to 
be described shortly, is our principal finding. The prospects 
for observing the transitions and the intermediate phase in an 
atomic gas will be discussed in the conclusion. 

Let us begin by discussing the phases in Fig. [T]in qualita- 
tive terms. The superfiuid transition of scalar bosons in two 
dimensions is of the KT type, mediated by the binding of vor- 
tices, suggesting that we consider the analogous defects of a 
polar condensate. In the polar state we may write <f> = ne , 
where n is a real unit vector and 8 is a phase variable, (and 
we have we set the density equal to unity). In this represen- 
tation taking n —> — n and 9 — s- 8 + tt maps the spinor to 
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FIG. 2: (Color online) A half vortex pair connected by a domain 
wall (shaded region) in the n vector (shown red) arising from easy 
axis anisotropy (the easy axis is horizontal). The blue arrows denote 
the configuration of the phase 9. 



itself. Thus the elementary vortex has circulation , or one 
half of the usual quantum of circulation, and coincides with a 
'disclination' in the vector n. 

The character of these point defects is dramatically altered 
by the inclusion of the Zeeman energy, which has the form 



H z = drtf [ P S z +qS 2 z ]cj>. 



(2) 



Here p and q describe the linear and quadratic effects (q > 0). 
We will be concerned with a system of zero total S z , so that 
the linear term has no effect: the case of non-zero S z will be 
discussed briefly at the end. In the (n, 6) representation the 
quadratic effect contributes an energy per particle of q(l — r? z ) , 
amounting to an easy-axis anisotropy for the n variable. 

At this point, it is convenient to introduce a simple lattice 
model which will be useful for our numerical simulations: 



H = H t + H z 



H 7 



'E'< 



H t = — ^4"l ' 4>j + C - C J = — 2f rii ■ iij cos (9i — 9j) 

(3) 
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with hopping parameter t. The corresponding continuum 
model for the spinor takes the form 



H 



d d r t(Vnf +t (V6>) 



qn. 
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(we take the lattice spacing equal to unity). Notice that the su- 
perfluid and magnetic degrees of freedom appear to decouple 
in this expression. The only coupling is global, in that half- 
vortex / disclinations are allowed. Thus when these defects 
are absent (or bound) the n degrees of freedom are described 
by the familiar Heisenberg model with anisotropy. 

In Eq. q appears as a 'mass' for deviations from the 
easy axis, meaning that such deviations are confined to 'do- 
main walls' of thickness oc q^ 1 / 2 that have an energy per unit 
length - or tension - oc q x l 2 . Furthermore, these domain walls 
can terminate on the half-vortices discussed above, see Fig. [2] 

Now we are in a position to understand the structure of 
the phase diagram. At q = 0, the only finite temperature 



transition is a KT transition at T = Tkt ~ tir/8 due to 
binding of the half-vortices, which have lower energy than 
integer vortices |7]|8). The distinctive features of this tran- 
sition will be discussed later. After the half-vortices become 
bound, the description of the n degrees of freedom coincides 
with that of the ordinary Heisenberg model, and no order ap- 
pears at finite temperature. For q nonzero but small the do- 
main walls connecting half-vortices have too little tension to 
affect the KT transition. Once the half-vortices have become 
bound, the domain walls are closed but otherwise fluctuating, 
and can disappear in an Ising transition at a lower tempera- 
ture. An alternative way to think about this second transition 
is in terms of the Heisenberg model, which has an Ising tran- 
sition in the presence of an easy-axis anisotropy. Finally, at 
large q, the n vector is pinned to the ±z direction and the 
model Eq. ([3) coincides with the usual XY model after shift- 
ing 9i — » 6i + nn Z: i/2. In this regime there is a single KT 
transition of integer vortices. In terms of the original spin 
states of the boson, only the S z = state is occupied, so that 
the behavior of a scalar condensate is recovered. An earlier 
investigation considered only q = and large q, so that the 
interesting interpolation between these two limits went unno- 
ticed Q. We note parenthetically that the case q < - harder 
to realize experimentally - was also studied recently (9). 

The region occupied by the intermediate PS phase in the 
model Eq. |3]l is very small: its presence could not be unam- 
biguously determined by Monte Carlo simulations on systems 
of up to 40 x 40 sites. This is likely due to the Ising transi- 
tion line being very steep near q = 0. Standard arguments for 
the scaling of the 'mass' (correlation length) with temperature 
ifTUl show that T ~ — l/log(g) for small q. We therefore 
study a generalized model with a 'pair hopping' term 
u ■ 



(5) 



(ij) 



Including H u in Eq. we see that a finite u only stiffens the 
phase variable, changing the coefficient of the 6 term from t 
to t + 2u. The half-vortex KT transition then occurs at the 
higher temperature T ss (t + 2u)ir/8, increasing the size of 
the PS phase. In the following we take u = 2t as then the PS 
phase is clearly visible even for moderate system sizes. Ex- 
perimentally a similar result could be achieved by increasing 
C2 until 2-body singlet bound states form at q = 0. Then n 
is disordered even at T = (enlarging the PS phase) until q 
is large enough to cause a quantum phase transition into the S 
phase. Photoassociation data suggest that the required condi- 
tion, a divergent singlet scattering length, has been achieved 
already via optical Feshbach resonance (Fig. 7 of Ref. ifTTI ). 

For q —> oo the model is equivalent to the Hamiltonian of 
the generalised XY model (see e.g.j2]3) 



E 

(ij) 



A cos(6>i - 9j) + (1 - A) cos(26». ( - 26 j 



(6) 



where conventionally A e [0,1]. This model also exhibits a 
PS phase. Our choice of parameters corresponds to the case 




FIG. 3: (Color online) The specific heat, C, for L — 16, showing 
the merging of the peaks associated with the two transitions as q 
increases. Inset: C for different system sizes at q = 1. 



A = 0.5 so that, according to the phase diagram in Ref. |4| 
the large q limit of our model with u = 2t still has a single 
transition (hence u itself does not produce a PS phase for q — > 
oo). As q is reduced this transition should split in two. 

We study the phase diagram of the model H = H t + H u + 
Hz via Monte Carlo simulations (using tools from the ALPS 
libraries [12]) on square systems of L x L sites with periodic 
boundary conditions. As there are three continuous param- 
eters per site a large number of sweeps of the lattice are re- 
quired to equilibrate and collect reliable data, even for small 
system sizes (> 10 6 for L = 8). We performed simula- 
tions for L = 8, 16, 24, 32 with some extra data collected for 
L = 40, 48 in special cases. To detect two separate transi- 
tions, we consider the specific heat capacity in addition to the 
Binder cumulants for both the spinor and the phase. 

A phase transition of the Ising type at T = T c should 
present itself as a sharp peak in the specific heat, C, for fi- 
nite size simulations. KT transitions are also accompanied by 
a peak in C, above Tkt, associated with the increase in en- 
tropy when vortices unbind. Fig. [3]shows that as q approaches 
a critical value, q ~ 4, the sharper lower temperature Ising 
peak and the broader higher temperature KT peak merge. 

Better quantitative information is given by the Binder cu- 
mulants ifTSl . The Binder cumulant for the spinor is 



Bl 



(7) 



where $ = J^. <j>i/N. An example plot is given in Fig. |4] We 
also calculate the cumulant for the z component: 



B{<f> z 



(8) 



These cumulants are sensitive to order and Ising-like order in 
(f>, respectively. On the other hand they are not sensitive to 



FIG. 4: (Color online) Intersection of the Binder Cumulants, B[ 
for different L at q — 1. Inset: enlarged region of intersection. 
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FIG. 5: (Color online) Left: Intersection of the Binder Cumulants 
B[<j>] for different L at q = 1. Right: The same data but plotted 
against (T — T C )L. Note that the data for L = 8 is not within the 
range of validity for finite size scaling. 



order (or quasi long-range order) in the phase alone. Instead, 
to look for order in exp(2i6>), we use the cumulant 

mm gtf*j (exppifr-flj+gfc-fl,)]) 
m ~ (E^exp»-^-)]» 2 ■ (> 

In the vicinity of a conventional, continuous transition at 
T c , and where finite size scaling holds, the Binder cumulant 
for a suitable variable can be written in the form 



Bl = B(TL"' 



(10) 



where B is a universal scaling function, and T = (T—T c ) /T c . 
From this we conclude that Binder cumulants for different L 
cross at T = 0, providing an accurate method for determining 
T c . For KT transitions, eq. ( 10 1 does not hold. However the 



crossings for different L still occur in a suitably narrow range 
lfl4ll (see Fig. [6]), allowing us to estimate Tkt- Fig. [TJshows 
T c and Tkt as found using the above cumulants. The values 
of T c provided by B[<j>] and B\(f> z \ agree within error over the 
full range of q investigated. As a further check on the nature 
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FIG. 6: (Color online) Intersection of the Binder Cumulants B[20] 
for different L at q — 1. Inset: enlarged region of intersection. 

of the PS to S transition we extract the critical exponent v by 
calculating dBr,/ dT at T c . We find that v is at least consistent 
with the Ising value v = 1, within error. An example of the 
resulting data collapse is shown in Fig. [5] 

Following Ref. Q we also examine the helicity modulus 
(also known as spin stiffness or superfluid density), T, de- 
fined as the change in free energy due to a twist in boundary 
conditions along some direction, x. A KT phase transition at 
Tkt is reflected in the helicity in the form of a jump propor- 
tional to Tkt (see e.g. ifTUl ). When the transition is driven by 
half-integer vortices this jump will be four times larger than 
for the integer case |fl5ll : 

ATi KT = AAT KT = (11) 

2 7T 

The observed position of the helicity jump as a function of 
q and T confirms that the transition temperature provided 
by B[26] represents Tkt with reasonable accuracy. As in 
Ref. Q we find that the N to PS transition is due to the 
presence of half vortices (Fig. [7]). In fact, setting u = 2t 
means that the transition is facilitated by half vortices even 
for q — » oo. This is consistent with eq. (|6]l for A = 0.5. 
In that case half-vortices still exist, though they occur in the 
phase alone and the line defects that join them have energy 
independent of q [4|. As m approaches zero the integer KT 
transition at large q should be recovered. 

We now return to the case of finite S z . As in the case of an 
antiferromagnet with an easy axis anisotropy, increasing S z 
leads to a spin-flop transition between a state with n aligned 
in the z -direction to one where it lies in the x — y plane. Such 
a transition is described by a bicritical point of the Heisenberg 
type, which in d = 2 must occur at T = H6] [17). At 
finite T, the low S z Ising and high S z , xy ordered, phases are 
separated by a normal region. The high S z phase resembles 
the q < case discussed in Ref. J9). 

In conclusion we have argued that polar condensates un- 
dergo separate KT and Ising transitions when subjected to the 



FIG. 7: (Color online) Helicity for q — 1. The upper and lower 
dashed lines are — and — respectively. The jump is clearly com- 
mensurate with AT = — . 



quadratic Zeeman effect. We have supported this finding with 
Monte Carlo simulations. Aside from the thermodynamic sig- 
natures discussed here, the PS phase should be visible in the 
correlation function of occupancies of different momentum 
states, as measured by the noise correlations in an image of 
the expanded gas: (<5n(ki)<5n(k2)) oc |ki + k2 | 4r?— 2 with 
v = T/(27tT) GO. 
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